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Abstract. - A recently introduced particle-based model for fluid dynamics with continuous ve- 

Slocities is generalized to model immiscible binary mixtures. Excluded volume interactions between 
the two components are modeled by stochastic multiparticle collisions which depend on the local 
'-^j . velocities and densities. Momentum and energy are conserved locally, and entropically driven phase 

separation occurs for high collision rates. An explicit expression for the equation of state is de- 
rived, and the concentration dependence of the bulk free energy is shown to be the same as that 
O ■ of the Widom-Rowlinson model. Analytic results for the phase diagram are in excellent agreement 

with simulation data. Results for the line tension obtained from the analysis of the capillary wave 
spectrum of a droplet agree with measurements based on the Laplace's equation. The introduction 
■ of "amphiphilic" dimers makes it possible to model the phase behavior and dynamics of ternary 

surfactant mixtures. 



o 

Introduction. — Hydrodynamic interactions and thermal fluctuations play a crucial 
role in a wide range of phenomena in soft matter physics and molecular and cellular biology. 
Because of the complexity of these systems, simulations have played an essential role in much 
of the research in these areas. In fact, the wide range of length and time scales in these 
problems places severe requirements on simulation protocol, and has lead to the development 
of several new coarse-grained, mesoscale simulation techniques such as lattice gas automata 
[1], the lattice Boltzmann method [2], dissipative particle dynamics [3,4], smoothed particle 
hydrodynamics [5] , and a newer approach variously called multi-particle collision dynamics 
or stochastic rotation dynamics (SRD) [6]. The basic motivation of all these approaches 
is to coarse-grain out irrelevant atomistic details while correctly incorporating the essential 
physics and conservation laws. 

SRD has several attractive features which have lead to its use in studies ranging from sed- 
imentation in colloidal suspensions [7] to the dynamic behavior of polymers in solution [8,9] 
and vesicles in flow [10]. In particular, it enables simulations in the microcanonical ensem- 
ble while fully incorporating both hydrodynamic interactions and thermal fluctuations; in 
addition, because SRD is a particle-based method, the coupling to colloidal particles, poly- 
mers, or other aggregates is straightforward, and the Brownian motion of these embedded 
objects is realized in a very natural way — through random collisions with the solvent parti- 
cles. Finally, the simplicity of the algorithm has made it possible to obtain accurate analytic 
expressions for the transport coefficients [11-13]. 
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Recently, it has been shown how to generalize the multi-particle collisions of the original 
SRD algorithm to model excluded volume effects, allowing for a more realistic modeling of 
dense gases and liquids with a nonidcal equation of state [14,15]. In this Letter we show that 
a similar approach can be used to model immiscible binary mixtures. The resulting model 
retains much of the simplicity of the original SRD algorithm, allowing for a detailed analysis 
of the transport coefficients and an explicit calculation of the equation of state and the bulk 
entropy and free energy densities. Since there is no potential energy, all interactions in the 
model are entropic, and the resulting bulk free energy density is the same as that of the 
Widom-Rowlinson model [16]. Theoretical predictions for the phase diagram are shown to 
be in good agreement with simulation data, and results for the line tension obtained from 
an analysis of the spectrum of capillary wave fluctuations agree with measurements based 
on the Laplace equation. 

Model. — In a binary mixture of A and B particles, phase separation can occur when 
there is an effective repulsion between A-B pairs. In the current model, this is achieved by 
introducing multi-particle collisions between A and B particles. The fluid is modeled by a 
large number TV of point-like particles of unit mass which move in continuous space with 
a continuous distribution of velocities. There are Na and Nb particles of type A and B, 
respectively In two dimensions, the system is coarse-grained into (L/a) 2 cells of a square 
lattice of linear dimension L and lattice constant a. The generalization to three dimensions 
is straightforward. 

To define the collisions, we introduce a second grid with sides of length 2a which groups 
four adjacent cells into one "supercell". As discussed in [17], Galilean invariance requires 
that the collisions occur in a randomly shifted grid. All particles are shifted by the same 
random vector with components in the interval [—a, a] before the collision step. Particles 
are then shifted back by the same amount after the collision. To initiate a collision, pairs of 
cells in every supercell are selected at random. As shown in Fig. 1, three different choices 
are possible: a) horizontal (with <j\ = x), b) vertical [cr-i = y), and c) diagonal collisions 
(with er 3 = (x + y)/ \pl and <T4 = {x — y)/ s/2) . For each pair of cells, two types of collisions 
are possible. As illustrated in Fig. 1, particles of type A in the first cell can undergo a 
collision with particles of type B in the second cell; vice versa, particles of type B in the 
first cell can undergo a collision with particles of type A in the second cell. There are no 
A-A or B-B collisions, so that there is an effective repulsion between A-B pairs. The rules 
and probabilities for these collisions are chosen in the same way as in the nonideal single- 
component fluid described in Refs. [14,18]. For example, consider the collision of A particles 
in the first cell with the B particles in the second. The mean particle velocity of A particles 
in the first cell is = (1/nu) v »i where the sum runs over all A particles, n%A, in 

the first cell. Similarly, = (1/tob) Y^f\ v « ^ s the mean velocity of B particles in the 
second cell. The projection of the difference of the mean velocities of the selected cell-pairs 
on (jj, Auab = &j • (ua - ub), is then used to determine the probability of collision. If 
Auab < 0, no collision will be performed. For positive Auab, a collision will occur with an 
acceptance probability 

PA(m Al m B , Auab) = max{l,imjm B Au A b Q(Au A b)} , (1) 

where is the unit step function and A is a parameter which allows us to tune the equation 
of state; in order to ensure thermodynamic consistency, it must be sufficiently small that 
AniAmB Auab ©(Au^b) < 1 for essentially all collisions [15]. When a collision occurs, the 
parallel component of the mean velocities of colliding particles in the two cells is exchanged, 
v\ (t + t) — u AB = — (vj (t) — u AB ), where u AB = (mAU A + mBU B )/(mA + mB) is the parallel 
component of the mean velocity of the colliding particles. The perpendicular component 
remains unchanged. It is easy to verify that these rules conserve momentum and energy in 
the cell pairs. The collision of B particles in the first cell with A particles in the second is 
handled in a similar fashion. 
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Fig. 1: Schematic binary collision rules. Momentum is exchanged in three ways: a) horizontally 
along <ti, b) vertically along 01, and c) diagonally along crz and <T4. w and Wd denote the proba- 
bilities of choosing collisions a), b) and c), respectively. For a cell pair defined by the vector cr, A 
particles in one cell collide with B particles in the other cell, and vice versa. 

Fig. 2: Non-ideal contribution to the pressure, P„, (in units of ksT/a 2 ) as a function of the tuning 
parameter A (in units of r/a). The bullets (•) and open circles (o) are results for ksT = 0.006 and 
ksT — 0.003, respectively. The lines are plots of Eq. (7). The inset shows the pressure difference 
across a droplet interface as a function of the inverse droplet radius, for A = 0. 60 at fcsT = 0.0005. 
The pressure is measured using the expression given in Eq. (8). The solid line is a plot of Eq. 
(13), using the line tension obtained from the analysis of capillary waves. Parameters: L/a — 64, 
M A = N A /(L/a) 2 = 5, M B = Nb /(L/a) 2 = 5, a = 1, and r = 1.0. 



Because there are no A-A and B-B collisions, additional SRD collisions at the cell level 
are incorporated in order to mix particle momenta. The order of A-B and SRD collision is 
random, i.e., the SRD collision is performed first with a probability of one half. If necessary, 
the viscosity can be tuned by not performing SRD collisions every time step. The results 
presented in this Letter were obtained using a SRD collision angle of 90°. 

Transport Coefficients. — The transport coefficients can be determined using the 
same Green-Kubo formalism as for the original SRD algorithm [11,17,19] and its extensions 
[14,15,18]. In particular, the discrete Green-Kubo relation 

= N¥T £ 9 a p(knr) , (2) 

B n=0 

with g a ff(k; t) = (k\o a \(0)\k\' o py (t)}, expresses the matrix of viscous transport coefficients, 
A a p (k) , in terms of a sum of time correlation functions of the reduced fluxes 

v ia (t)t ■ A^(t) + Av ja (t)k ■ [ACjit) - z s jt (t + r)/2] - ^«?(t)j , 

(3) 

where (t) = ^ (t + r) - ^ (t), Av jx (t) = v jx (t + r)- v jx (t) and (t) - ^ (t + r) - 
£j (t + t) . t is the time step in the simulation and the prime on the sum in Eq. (2) indicates 
that the n = term has the relative weight 1/2. £j(t) is the cell coordinate of particle j 
at time t = tit, while is it's cell coordinate in the (stochastically) shifted frame, ^{t) 
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indexes pairs of cells which participate in a collision event at time r; the second subscript, 
I, is the index of the collision vectors rr; shown in Fig. 1. For example, for collisions 
characterized by a±, Zj lx = 1 if £!j x in Eq. (3) is one of the two cells on the left of a supercell 
and Zj lx = — 1 if £j X is on the right hand side of a supercell; all other components of are 
zero. In general [18], the components of are either 0, 1, or — 1. 

Assuming only cubic symmetry, the most general form of A Q!| a(k) in two dimensions 
is [11] 

A Q/3 (k) = vi5 a p + v 2 (5 a p - k a kf^J + ~/k a k f3 + e k-yk p I a p 7P , (4) 

where / is the rank four unit tensor, vi is a new viscous transport coefficient associated 
with the non-symmetric part of the stress tensor, 7 is the bulk viscosity and e is a viscosity 
coefficient which is nonzero only if there is cubic anisotropy. In a simple liquid, e = 
(because of invariance with respect to infinitesimal rotations), v = z^i, and ^2 = (because 
the stress tensor is symmetric in d a vp). 

Since both particle streaming and the collisions contribute to momentum transport, 
there are two — kinetic and collisional — contributions to the transport coefficients. The ki- 
netic contribution dominates at large mean free paths, A, the collisional for A/a -C 1. For 
the original SRD algorithm, the kinetic contribution is isotropic, so that there is only one 
viscosity, v\ the kinetic contribution to the bulk viscosity is also zero, as in a real ideal 
gas [13]. However, because SRD collisions do not, in general, conserve angular momentum, 
the microscopic stress tensor is not symmetric and there is a collisional contribution [11] 
to vi- It should be noted, however, that the SRD algorithm can be modified to conserve 
angular momentum in two dimensions by introducing a unique, configuration dependent 
collision angle in each cell [20]. 

In Rcfs. [14,18] it was argued (for the nonideal model) that the probabilities of horizontal 
and vertical (w) and diagonal (wd) collisions should be chosen so that the relaxation rate of 
the second moments of the velocity distribution function all decay at the same rate. This 
lead to the requirement that w = 1/4 and Wd — 1/2. Here we show that the same result 
follows from the requirement that the kinetic contribution to the viscous stress tensor is 
symmetric, so that there is only one viscosity, and e = 7 — = 0. To see this, note first that 
for a = (3 = 1 with k = y and k = x, Eqs. (2)-(4) give 7 - v 2 = -e/2 = [An(y) - Au(f)] 
and v\ + V2 = An(y). For large mean free paths, the term in parentheses in Eq. (3) 
reduces to Vj a (t)\z- Vj(t) — (rfc Q /2)w|(i) in two dimensions; using this result and making the 
assumption of molecular chaos, the sum in Eq. (2) reduces to a geometric series in .gn(k; r). 
The calculation of these quantities requires that the contributions from the horizontal, g^[, 
vertical, g\ 1: and diagonal collisions, g^[, be calculated individually and then summed in 
the form 

g 11 (Ur)=w[gf 1 (k;T)+gY 1 (kr)]+w d g^,T) . (5) 

If <7ii(j/;t) ^ gii(x;r), there is a cubic anisotropy. Ignoring fluctuations in the number of 
particles per cell, we find that 511(2/; t) = on(x;r) only if w = 1/4 and Wd = 1/2. In this 
case, the only non-zero viscosity is 

v ^K 11 ^) = T J^(^^^[M A M B {M A + M B )]-^-l^ , (6) 

where Ma = ^/(L/a) 2 and Mb = Nb/{L/ci) 2 . In deep quenches, the density of the 
minority phase is very small, and the nonideal contribution to the viscosity approaches zero, 
i.e. g(r) ~ 1 — 0(p A )\ in this case, the SRD collisions provide the dominant contribution to 
the viscosity. 

Free Energy. — An analytic expression for the equation of state of this model can be 
derived by calculating the momentum transfer across a fixed surface, in much the same way 
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Fig. 3 Fig. 4 

Fig. 3: Phase diagram of a 50% A - 50% B mixture. There is phase separation for pY > 2. The 
inset shows a configuration 50, 000 time steps after a quench along pab = to pT — 3.62. The dark 
(blue) and light (white) colored spheres are A and B particles, respectively. Parameters: L/a — 64, 
M A = M B = 5, k B T = 0.0004, r = 1 and a = 1. 

Fig. 4: Dimensionless radial fluctuations, (|rjb | 2 ), as a function of the mode number k for yl = 0.45 (•) 
and A = 0.60 (o), with ksT = 0.0004. The average droplet radii are ro = 11.95 a and ro = 15.21 a, 
respectively. The solid lines are fits to Eq. (12). The inset shows a typical droplet configuration 
for A = 0.60 and T = 0.0004. Parameters: L/a = 64, M A = 2, Mb = 8, a = 1 and t = 1. 



as was done for the nonideal model in Ref. [14]. The resulting expression for the nonideal 
contribution to pressure is 

P n =(w+^j AMaMb 1 ^- = Tk B T PAPB: (7) 

where pa and p B are the densities of A and B and r = (w + Wd/V2)a 3 A/r. In simulations, 
the total pressure can be measured by taking the ensemble average of the diagonal compo- 
nents of the microscopic stress tensor. In this way, the pressure can be measured locally, at 
the cell level. In particular, the pressure in a region consisting of N c cells is 

Av jx z? lx /2]\ , (8) 

where the second sum runs over the particles in cell c. The first term in Eq. (8) is the 
ideal gas contribution to the pressure; the second term comes from the momentum transfer 
between cells involved in the collision indexed by z s . The results of measurements of the 
non-ideal contribution to the pressure obtained using Eq. (8) are shown in Fig. 2 for 
k B T = 0.006 (•) and k B T = 0.003 (o). The lines are the theoretical predictions of Eq. (7). 
For small values of A, the agreement between theory and simulation is excellent; deviations 
for larger A are an indication that the model is no longer thermodynamically consistent [15]. 

Equation (7) can be used to determine the entropy density, s. The ideal gas contribution 
to s has the form [21] 

Sideai = ptp(T) -k B [p A In p A + p B In p B ] , (9) 

where p — pa + p B - Since (f(T) is independent of pa and p B , this term does not play a 
role in the current discussion. The nonideal contribution to the entropy density, s n , can be 



1 Nc 

\c— 1 i£c 
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obtained from Eq. (7) using the thermodynamic relation [21] 

P n /T = s n - p A ds n /dp A - p B ds n /dp B . (10) 

The result is s n = —k B Tp A p B , so that the total configurational contribution to the entropy 
density is 

s = -k B [pA hi Pa + Pb hips + Tp A p B ] . (11) 

Since there is no configurational contribution to the internal energy in this model, the 
mean field phase diagram can be determined by maximizing the entropy at fixed density 
p. The resulting demixing phase diagram as a function of p AB = (pa ~ Pb)/ P 1S given 
by the solid line in Fig. 3. The critical point is located at pab = 0, (pT)* = 2. For 
pT < 2, the order parameter p AB = 0; for pT > 2, there is phase separation into coexisting 
A and B-rich phases. Simulation results for p AB obtained from density histograms are 
shown as bullets (•). The dashed line is a plot of the leading singular behavior, p AB = 
y/3(pT — 2)/2, of the order parameter at the critical point. As can be seen, the agreement 
between the mean field predictions and simulation are very good except close to the critical 
point, where the histogram method of determining the coexisting densities is unreliable and 
critical fluctuations could influence the shape of the coexistence curve. 

Line Tension. — A configuration after 50,000 time steps following a quench to point 
Pab = 0, pT = 3.62 of the phase diagram is shown in the inset to Fig. 3, and a snapshot 
of a fluctuating droplet at p AB = —0.6, pT = 3.62 is shown in the inset to Fig. 4. The 
amplitude of the capillary wave fluctuations of a droplet is determined by the line tension, 
a. In particular, for a droplet in an incompressible fluid, the mean square amplitude of 
fluctuations of the Fourier components of the (dimensionless) droplet radius, (|rfc| 2 ), is 
related to the line tension by the dispersion relation [22] 

(hi 2 ) - ^ (J^) , (12) 

where ro is the mean radius of the droplet. Fig. 4 contains a plot of (|?"fc| 2 ) as a function 
of mode number k for A = 0.45 (pT = 3.62) and A = 0.60 (pT = 2.72). Fits to the data 
yield aa/k B T ~ 2.9 for pT = 3.62 and aa/k B T ~ 1.1 for pT = 2.72. Mechanical equilibrium 
requires that the pressure difference across the interface of a droplet satisfies the Laplace 
equation, 

Ap = p in -p out = a/r . (13) 

Using Eq. (8), we have determined Ap as a function of the droplet radius for A = 0.60 and 
k B T = 0.0005. The results are plotted in the inset to Fig. 2, where it can be seen that the 
Laplace relation is satisfied for the correct values of the line tension. 

The model therefore displays the correct thermodynamic behavior and interfacial fluc- 
tuations. It can also be extended to model amphiphilic mixtures by introducing dimcrs 
consisting of tethered A and B particles [22]. If the A and B components of the dimcrs 
participate in the same collisions as the solvent, they behave like amphiphilic molecules in 
oil- water mixtures. The resulting model displays a rich phase behavior as a function of pT 
and the number of dimers, N4. We have observed both droplets and micelles, as shown in 
Fig. 5, and a bicontinuous phase, as illustrated in Fig. 6. The coarse-grained nature of the 
algorithm enables the study of large time scales with a feasible computational effort. 

* * * 

Support from the National Science Foundation under Grant No. DMR-0513393 and ND 
EPSCoR through NSF grant EPS-0132289 are gratefully acknowledged. 
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Fig. 5 Fig. 6 



Fig. 5: A snapshot showing a droplet in a mixture with Na = 8, 192, Nb = 32, 768 and N& = 1, 500 
dimers after 10 time steps. The initial configuration is a droplet with a homogeneous distribution 
of dimers. The dark (blue) and light (white) colored spheres indicate the A and B particles, 
respectively. For clarity, A particles in the bulk are smaller, and B particles in the bulk are not 
shown. Parameters: L/a = 64, M A = 2, Mb = 8, A = 1.8, k B T = 0.0001, r = 1 and a = 1. 

Fig. 6: Typical configuration showing the bicontinuous phase for Na = Nb = 20, 480 and Nd = 
3, 000. Parameters: L/a = 64, M A = 5, M B = 5, A = 1.8, k B T = 0.0001, r = 1 and a = 1. 
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